Estimating the time of arrival of a seismic wave

ABSTRACT

Data processing means are described for calculating an estimated time of arrival of a seismic or microseismic P or S wave at a sensor station, based on a P to S wave velocity ratio, a calculated estimated time of origin of the event and, where the estimated arrival time of a P wave is to be calculated, a picked arrival time of the S wave at the sensor station or, where the estimated arrival time of a S wave is to be calculated, a picked arrival time of the P wave at the sensor station.

The present invention relates to a method of calculating an estimatedtime of arrival of a seismic or microseismic P or S wave at a sensorstation, data processing means for carrying out such a method, and adata carrier comprising computer readable program means for carrying outsuch a method when run on a computer.

A typical seismic or microseismic event produces both P-wave and S-wavesignals. These distinct signals travel at different velocities throughthe earth. For a wide range of crustal rock types the ratio between theP-wave and S-wave velocity is well known and does not vary to any greatextent. Indeed in many instances it is not possible a priori to obtainan-accurate S-wave velocity and so the S-wave velocity structure for aregion of interest is often derived by assuming a constant ratio withthe P-wave velocity structure. Thus we find that, in practice, the twoseismic velocities structures within a volume of interest will show anapproximately constant ratio to each other, even where the value of bothtypes of velocity changes.

The normal method of estimating the location and time (in terms of thecoordinates x, y, z and time zero (t0)) at which a seismic ormicroseismic event occurs involves positioning a number of sensorstations, comprising geophones or the like, at various locations in acatchment area, and measuring the arrival times at each of said stationsof the P and S waves generated. Once the time of arrival of each wave ateach station has been picked (i.e. decided upon based on the seismicdata recorded at the station concerned), a guess is made at the locationand time of the seismic or microseismic event, and the known or assumedP and S wave velocities are used to calculate the times at which the Pand S waves would reach each station, assuming that this guess werecorrect. These calculated arrival times, based on the guessed locationand time of origin, are then compared with the picked arrival times(which, as noted above, are based on the seismic data actually recorded)and the location and time of origin adjusted, via a mathematicalprocesses known as the least squares process, such that the finalestimated location and time of origin of the event best fits, in a leastsquares sense, with the picked arrival times of the waves at thestations.

While the above method is generally accepted to be the best method ofevent location currently available, the reliability and accuracy of thelocation of a given seismic or microseismic event, recorded on a givennetwork, is heavily dependent upon the accuracy and reliability of thearrival time picking. In some cases, the quality of the seismic data maymake it difficult to identify a distinct arrival time at all. If apicked arrival time is badly incorrect, or worse still a wave typewrongly identified (for example where a P-wave picked as an S-wave or aP to S conversion picked as an S-wave), the final estimated location andtime of origin of the event is likely to differ significantly from thetrue location and time. Furthermore, the least squares location processmeans that the effect of a single bad timing is spread out and manifestsitself in all the channel residuals and the overall misfit of thelocation, making the identification of incorrect pick less thanstraightforward.

According to one aspect of the present invention, data processing meansfor calculating an estimated time of arrival of a seismic ormicroseismic P or S wave at a sensor station are provided, comprising adata processor adapted to:

-   -   a) calculate an estimated time of origin for the seismic or        microseismic event generating the P and S waves, based on a P to        S wave velocity ratio and picked arrival times of the P and S        waves at a sensor station other than the one for which the        estimated time of arrival of the P or S wave is to be        calculated; and    -   b) calculate the estimated time of arrival of the P or S wave,        based on a P to S wave velocity ratio, the estimated time of        origin of the event and, where the estimated arrival time of a P        wave is to be calculated, a picked arrival time of the S wave at        the sensor station for which the estimated arrival time of the P        wave is being calculated or, where the estimated arrival time of        a S wave is to be calculated, a picked arrival time of the P        wave at the sensor station for which the estimated arrival time        of the S wave is to be calculated.

The data processor may be adapted to calculate estimated arrival timesfor both the P and S waves at a sensor station.

Preferably, the data processor is adapted to calculate a plurality ofestimated times of arrival of the P and/or S wave at a sensor station,based on a plurality of estimated times of origin for the microseismicevent calculated from the picked arrival times of the P and S waves at aplurality of sensor stations other than the one at which the estimatedtimes of arrival are to be calculated. The data processor may be furtheradapted to display the picked arrival times and estimated arrival timesin relation to each other such that the clustering pattern of thearrival times can be analysed. Where this is the case, the dataprocessor may also be adapted to display information regarding thecalculation of any particular estimated arrival time in response to theselection of said estimated arrival time by an user.

Preferably, the data processor is adapted to calculate one or moreestimated times of arrival for the P and/or S waves at each sensorstation in a network of sensor stations, wherein for each sensor stationthe necessary estimated time or times of origin are calculated from thepicked arrival times of the P and S waves at one or more of the otherstations in said network.

In some cases, the data processor may require the user to pick arrivaltimes for the P and S waves at the various sensor stations by studyingthe seismic data recorded at each station, and to then provide the dataprocessor with the picked arrival times. Alternatively, the dataprocessor may be adapted to receive seismic data from the sensorstations and to pick arrival times for the P and S wave at each sensorstation based on said seismic data.

Where a number of possible arrival times for a wave at a sensor stationcould be picked (either by the user, or by the data processor itself),the data processor may adapted to compare said possible arrival timeswith any estimates calculated for the arrival time of said wave at saidstation in order to determine which of the possible picked arrival timesare more likely to correspond to the true arrival time of said wave atsaid sation. The data processor may be further adapted to select ormodify one of said possible arrival times in order to arrive at a finalpicked arrival time that, based on the above determination, seems mostlike to correspond to the true arrival time of said wave. Alternatively,the data processor may be adapted to indicate which of the possiblearrival times should be selected or modified in order to arrive at afinal picked arrival time that seems, based on the above determination,to be most likely to correspond with the true arrival time.

According to another aspect of the present invention, a data carriercomprising computer readable program means for adapting a computer tofunction as the data processing means according to the above aspect ofthe invention is provided.

According to a further aspect of the present invention, a method ofcalculating an estimated time of arrival of a seismic or microseismic Por S wave at a sensor station is provided, said method comprising thesteps of:

-   -   a) calculating an estimated time of origin for the seismic or        microseismic event generating the P and S waves, based on a P to        S wave velocity ratio and picked arrival times of the P and S        waves at a sensor station other than the one for which the        estimated time of arrival of the P or S wave is to be        calculated; and    -   b) calculating the estimated time of arrival of the P or S wave,        based on a P to S wave velocity ratio, the estimated time of        origin of the event and, where the estimated arrival time of a P        wave is to be calculated, a picked arrival time of the S wave at        the sensor station for which the estimated arrival time of the P        wave is being calculated or, where the estimated arrival time of        a S wave is to be calculated, a picked arrival time of the P        wave at the sensor station for which the estimated arrival time        of the S wave is to be calculated.

Thus the present invention provides means for obtaining one or moreestimates as to the arrival time of the P and/or S waves at one or moreof the sensing stations, without the location of the seismic ormicroseismic event first having to be located. The estimated arrivaltimes thus produced can be used in a number of ways. Where the seismicdata recorded at a station cannot be interpreted with sufficientreliability in order to make a pick, an estimate of the arrival time canbe used to replace the pick that is not being made. Where a number ofsensing stations are present, a spread of estimated arrival times can becalculated, in order to provide the user with sufficient information toidentify possible bad picks before an attempt at event localisation ismade. Such spreads of estimates can also be used for analysis of the Pto S wave velocity structure and ratio, as part of an auto-locationprocess and/or for estimating picking uncertainty.

A specific embodiment of the invention will now be described, solely byway of example.

According to the present embodiment, the data processing means comprisea computer running software for calculating the necessary estimatedtimes of origin and arrival (though this could equally be achieved viasuitable electronic hardware) by applying certain algorithms, thederivation and application of which are explained below.

If one assumes that the network of sensor stations (each of which couldbe a seismic detector as described in GB-A-2 275 337 for example) aredeployed in an arbitrarily complex isotropic earth with a constant ornear to constant P- to S-wave ratio, then the distance (D) of amicroseismic event from a station can be expressed as:D=(Ts−To)Vsand/orD=(Tp−To)Vpwhere Vs is the S-wave velocity, Vp is the P-wave velocity, Tp is thearrival time of the P-wave at the station, Ts is the arrival time of theS-wave at the station, and To is the time of origin of the microseismicevent generating said P- and S-waves.

Multiplying through the two equations gives:TsVs−ToVs=Dand/orTpVp−ToVp=D

Eliminating D, the source to the distance, we get one equation:TsVs−ToVs=TpVp−ToVp

Collecting all terms including To, the origin time, gives:ToVp−ToVs=TpVp−TsVs∴ToVpVs−To=TpVp/Vs−Ts∴To(Vp/Vs−1)=TpVp/Vs−TsVp/Vs represents the P-wave to S-wave velocity ratio, which we willdesignate as R. Substituting R for Vp/Vs gives:To (R−1)=TpR−Ts∴To=(TpR−Ts)/(R−1)

We now have an estimate of the origin time of the earthquake in terms ofthe arrival times (Tp and Ts) of the P- and S-waves at the station andthe ratio R between the P-wave to S-wave velocity field ratio.

This equation has R in it twice and in order to re-arrange it so that Ronly appears once we can add a +Tp and a −Tp term to the top line of theequation, which gives:

${To} = \frac{{TpR} + {Tp} - {Tp} - {Ts}}{\left( {R - 1} \right)}$

Re-arranging we get:

${To} = \frac{{{Tp}\left( {R - 1} \right)} + {Tp} - {Ts}}{\left( {R - 1} \right)}$

Separating terms gives:

${To} = {\frac{{Tp}\left( {R - 1} \right)}{\left( {R - 1} \right)} + \frac{Tp}{\left( {R - 1} \right)} - \frac{Ts}{\left( {R - 1} \right)}}$

Now we can cancel the (R−1) term in the first part of the right-handside of the equation to give:

${To} = {{Tp} + \frac{Tp}{\left( {R - 1} \right)} - \frac{Ts}{\left( {R - 1} \right)}}$

To finally give:To=Tp−(Ts−Tp)/(R−1)

Thus, where the arrival times of the P- and S-waves at a station can bepicked with a reasonable degree of certainty and given a constantvelocity-ratio earth model we can, from that station, get an estimate ofthe origin time of the event based solely on the picked arrival times ofthe P- and S-waves and the assumed P- and S-wave velocity field ratio.Note that this equations does not require knowledge of the actuallocation of the microseismic event or the velocity fields for the P- andS-waves.

Where the station for which the P- and S-wave arrivals have been pickedis designated as station 1, we can add a 1 subscript to show that theseterms apply to station 1 to get:To ₁ =Tp ₁−(Ts ₁ −Tp ₁)/(R−1)  (1)where To₁ is the origin time estimated by using timings from station 1.

Likewise, if the arrival times were picked at a different station,designated station 2, we can add a 2 subscript to get:To ₂ =Tp ₂−(Ts ₂ −Tp ₂)/(R−1)  (2)

Note that, given the assumption of an isotropic earth, the velocityfield ratio remains constant for both stations.

However, the times of origin calculated from stations 1 and 2 should beidentical, since the same event is generating the waves recorded at bothstations. Thus, in circumstances where an S-wave pick at station 2 hasnot been possible, or the accuracy of the pick seems uncertain, suchthat an estimate for the S-wave arrival time at station 2 is desired,this can be obtained by substituting To₂ for To₁ and Ts_(e2) for Ts₂ inequation (2) to get:

   To₁ = Tp₂ − (Ts_(e2) − Tp₂)/(R − 1) $\begin{matrix}{{\therefore\mspace{11mu}{{To}_{1}\left( {R - 1} \right)}} = {{{Tp}_{2}\left( {R - 1} \right)} - \left( {{Ts}_{e2} - {Tp}_{2}} \right)}} \\{= {{{Tp}_{2}\left( {R - 1} \right)} - {Ts}_{e2} - {Tp}_{2}}}\end{matrix}$ ${\begin{matrix}{{\therefore\mspace{11mu}{Ts}_{e2}} = {{{Tp}_{2}\left( {R - 1} \right)} - {{To}_{1}\left( {R - 1} \right)} + {Tp}_{2}}} \\{= {{{Tp}_{2}R} - {Tp}_{2} - {{To}_{1}R} - {To}_{1} + {Tp}_{2}}} \\{= {{{Tp}_{2}R} - {{To}_{1}R} - {To}_{1}}} \\{= {{\left( {{Tp}_{2} - {To}_{1}} \right)R} - {To}_{1}}}\end{matrix}\;\therefore\mspace{11mu}{Ts}_{e2}} = {{R\left( {{Tp}_{2} - {To}_{1}} \right)} - {To}_{1}}$

Where Ts_(e2) represents an estimate for the S-wave arrival time atstation 2. Note that the calculation of the estimated time of arrival ofthe S-wave at station 2 requires only the P and S picks from station 1(in order to calculate To₁), the P-wave pick at sation 2 and the assumedconstant R ratio.

Likewise, where an estimate for the P-wave arrival time at station 2 isdesired, this can be obtained by substituting To₂ for To₁ and Tp_(e2)for Tp₂ in equation (2) to get:

   To₁ = Tp_(e2) − (Ts₂ − Tp_(e2))/(R − 1) $\begin{matrix}{{\therefore\mspace{11mu}{{To}_{1}\left( {R - 1} \right)}} = {{{Tp}_{e2}\left( {R - 1} \right)} - \left( {{Ts}_{2} - {Tp}_{e2}} \right)}} \\{= {{{Tp}_{e2}\left( {R - 1} \right)} - {Ts}_{2} - {Tp}_{e2}}}\end{matrix}$ $\begin{matrix}{{\therefore\mspace{11mu}{{{To}_{1}\left( {R - 1} \right)} + {Ts}_{2}}} = {{{Tp}_{e2}\left( {R - 1} \right)} + {Tp}_{e2}}} \\{= {{{Tp}_{e2}R} - {Tp}_{e2} + {Tp}_{e2}}} \\{= {{Tp}_{e2}R}}\end{matrix}$ ${\begin{matrix}{{\therefore\mspace{14mu}{Tp}_{e2}} = {\left( {{{To}_{1}\left( {R - 1} \right)} + {Ts}_{2}} \right)/R}} \\{= {\left( {{{To}_{1}R} - {To}_{1} + {Ts}_{2}} \right)/R}} \\{= {{{To}_{1}{R/R}} + {\left( {{Ts}_{2} - {To}_{1}} \right)/R}}}\end{matrix}\;\therefore\mspace{11mu}{Tp}_{e2}} = {{To}_{1} + {\left( {{Ts}_{2} - {To}_{1}} \right)/R}}$where Tp_(e2) is the estimated P-wave arrival time for station 2.

It follows, that where P-wave and S-wave picks for station 2 have beenmade, the process described above can be reversed in order to obtainestimates for the picks, both P-wave and S-wave, on station 1. Indeed,using the above process, every station for which a P-wave and an S-wavepick has been made allows an estimated time of arrival to be calculatedfor the P- and/or S-waves at another station in the network (assuming anarrival time for the P- or S-wave not being estimated has been pickedfor said station).

Thus, for a network of six stations, if all the P-wave and S-wave picksare made at each station then five estimated times of arrival for eachwave at each station can also be calculated. If a pick cannot be madefor the arrival time of one wave at one of the stations (due to thearrival time of the wave not being identifiable from the seismic data)then one of the five estimated arrival times for said wave at saidstation (or an average thereof) can be used for the purpose ofsubsequent localisation of the microseismic event. Likewise, if all thepicks have been made, but one of the picks is significantly differentfrom the true arrival time (again due to the difficult of identifyingthe true arrival time from the data recorded), then the estimatesproduced from the bad pick will differ significantly from the otherestimates which do not include the bad pick, allowing the bad pick to beidentified and corrected (or indeed replaced by an estimate) prior tothe localisation process. This is in contrast to attempting to examinethe channel residuals produced by the localisation process in order toidentify a bad pick, since each channel residue is influenced by all thepicks used in the location as well as the P to S velocity structureratio and three timings.

The estimated arrival times that have been generated can also be used tohelp improve the consistency of arrival time picking. If the assumptionof a constant P to S velocity ratio is true then all the estimatedarrival times for a wave should form a tight cluster around the pickedarrival time for a wave. If such a clustering pattern is not observed,there are three basic reasons why this may be the case, which are asfollows:

-   -   1) the assumption of a constant p to s ratio is does not hold        true for the area concerned;    -   2) the p to s ratio used is incorrect or    -   3) one or more of the picked arrival times is incorrect to a        significant degree.

Reasons (1) and (2) are related to the velocity structure. If the reasonis reason (1) then this can be easily remedied by changing the P or Swave velocity until good clustering is observed. If the reason is reason(2) then the P to S ratio can be adjusted until the optimum ratio isfound that gives the best clustering around a group of many timings frommany events. Thus the process of estimating arrival times describedherein could also be applied in order to optimise the velocitystructure. Typically it is the S-wave velocity that is most poorlyresolved and so in most cases it will be this velocity that needs to bealtered in order to obtain the optimum velocity ratio and structure. Itis possible to distinguish between reasons (1) and (2) in that wherereason (2) applies (ie. where the ratio is constant but inappropriate)good clustering together of the estimated picks will occur, but no thearound the picked arrival times. Where reason (1) applies the effectwill be more subtle and the estimated picks will not cluster to the sameextent. In practice this may be a more difficult situation to identify,but in most cases the variations in the P to S ratio within a particularregion are likely to be small, such that this will at most be a secondorder effect. Once many events have been located, it would be possibleto investigate the relationship between location and P to S ratio basedon the data obtained.

In most cases reason (3) will be the cause of the desired clusteringpattern not being observed. In order to assist the user in identifyingwhich stations picked arrival times are giving the bad estimate, thesoftware may further allow for the user to select any of the displayedestimates using the computer mouse or curser keys, in response to whichthe station that gave rise to the estimate is highlighted. In this waythe program implement on the computer guides the user to the stationwith the “worst” picks. This process can furthermore be repeated untilthe clustering of the estimated and measured picks appears satisfactoryto the user.

The program running on the computer may also be designed to use theestimated picks as part of an auto-location process. As has beendescribed above, a “good” pick can be identified by the fact that it hasa cluster of estimated picks around it. During the process ofauto-timing it is common for the computer to initially select a numberof potential picks for the arrival time of a wave at a station. Byapplying the above process of generating estimates for each pick andanalysing the clustering of the estimates around each prospectiveauto-picks, the program can be designed to eliminate those prospectiveauto-picks that are less likely to correspond to the true arrival time.The program can also use the degree of clustering to decide if an eventhas been picked well enough to warrant storing and use for a subsequentlocalisation process.

While estimating picking uncertainty or accuracy can be a difficult andsubjective process, it is also possible to use the estimated arrivaltimes in order to provide an indication of picking uncertainty. Eachestimated arrive time includes two other arrival times. For example, theS-wave arrive time on station n involves the following picks:

Tp₁ Ts₁ Tp_(n) Tp₂ Ts₂ Tp_(n) Tp₃ Ts₃ Tp_(n) Tp₄ Ts₄ Tp_(n) Tp₅ Ts₅Tp_(n)

In which Tp₁ indicates the picked arrival time for the P-wave at station1, Ts₁ indicates the picked arrival time for the S-wave at station 1,Tp_(n) indicates the picked arrival time of the P-wave at station n, andso on. From the above it will be apparent that Tp_(n) is a constant inall the estimates. Thus the scatter in the estimates of the Ts_(n) willbe a function of the picking uncertainty in the P- and S-wave pickingfor each channel.

Some of the above steps or components of the inventive method orapparatus are summarized in the block diagram of FIG. 1

1. Data processing apparatus for calculating an estimated time ofarrival of a seismic or microseismic P or S wave at a sensor station,the seismic or microseismic P or S wave being generated by a seismic ora microseismic event, comprising a data processor that: a) calculates anestimated time of origin for the seismic or the microseismic eventgenerating the P and S waves, based on a P to S wave velocity ratio andpicked arrival times of the P and S waves at a sensor station other thanthe one for which the estimated time of arrival of the P or S wave is tobe calculated; and b) calculates the estimated time of arrival of the Por S wave without knowledge of a location of the seismic or microseismicevent based on a P to S wave velocity ratio, the estimated time oforigin of the event and, where the estimated arrival time of a P wave isto be calculated, a picked arrival time of the S wave at the sensorstation for which the estimated arrival time of the P wave is beingcalculated or, where the estimated arrival time of a S wave is to becalculated, a picked arrival time of the P wave at the sensor stationfor which the estimated arrival time of the S wave is to be calculated.2. The apparatus according to claim 1, wherein said data processorcalculates estimated arrival times for both the P and S waves at asensor station.
 3. The apparatus according to claim 1, wherein said dataprocessor is configured to calculate a plurality of estimated times ofarrival of the P and/or S wave at a sensor station, based on a pluralityof estimated times of origin for the seismic or the microseismic eventcalculated from the picked arrival times of the P and S waves at aplurality of sensor stations other than the one at which the estimatedtimes of arrival are to be calculated.
 4. The apparatus according toclaim 3, wherein the data processor is further configured to display thepicked arrival times and estimated arrival times in relation to eachother such that the clustering pattern of the arrival times can beanalysed.
 5. The apparatus according to claim 4, wherein the dataprocessor comprises a display for displaying information regarding thecalculation of any particular estimated arrival time in response to theselection of said estimated arrival time by a user.
 6. The apparatusaccording to claim 1, wherein said data processor calculates one or moreestimated times of arrival for the P and/or S waves at each sensorstation in a network of sensor stations, wherein for each sensor stationthe necessary estimated time or times of origin are calculated from thepicked arrival times of the P and S waves at one or more of the otherstations in said network.
 7. The apparatus according to claim 1, whereinthe data processor receives seismic data from the sensor stations andpicks arrival times for the P and S wave at each sensor station based onsaid received seismic data.
 8. The apparatus according to claim 1,wherein, with a number of possible arrival times for a wave at a sensorstation, the data processor processes said possible arrival times withany estimates calculated for the arrival time of said wave at saidstation in order to determine an arrival time of said wave at saidstation.
 9. The apparatus according to claim 8, wherein the dataprocessor is configured to select or modify one of said possible arrivaltimes in order to arrive at a final picked arrival time based on saiddetermination.
 10. The apparatus according to claim 8, wherein the dataprocessor processes which of the possible arrival times should beselected or modified in order to arrive at a final picked arrival timebased on said determination.
 11. A method of calculating an estimatedtime of arrival of a seismic or microseismic P or S wave at a sensorstation, the seismic or microseismic P or S wave being generated by aseismic or a microseismic event, said method comprising the steps of: a)calculating an estimated time of origin for the seismic or themicroseismic event generating the P and S waves, based on a P to S wavevelocity ratio and picked arrival times of the P and S waves at a sensorstation other than the one for which the estimated time of arrival ofthe P or S wave is to be calculated; and b) calculating the estimatedtime of arrival of the P or S wave, based on a P to S wave velocityratio, the estimated time of origin of the event and, where theestimated arrival time of a P wave is to be calculated, a picked arrivaltime of the S wave at the sensor station for which the estimated arrivaltime of the P wave is being calculated or, where the estimated arrivaltime of a S wave is to be calculated, a picked arrival time of the Pwave at the sensor station for which the estimated arrival time of the Swave is to be calculated, wherein the estimated time of arrival iscalculated without knowledge of a location of the seismic ormicroseismic event.